##########################################################################################

library(data.table)
library(optparse)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/subcell/cluster17/cluster17.combineRNA.motif_peak2gene.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/motif_position"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
a <- load(comine_data_file)
## atac_proj

###########################################################################################
## 细胞类型
cell_type <- unique(atac_proj@cellColData$cell_type)
cluster <- unique(atac_proj@cellColData$seurat_clusters)

if( length(which(cluster == c(12,13))) == 2 ){
    cluster <- 17
}

###########################################################################################
## peak
GSM_se <- getMatrixFromProject(atac_proj, useMatrix="PeakMatrix")
GSM_mat <- assays(GSM_se)$PeakMatrix
tmp_name <- data.frame(GSM_se@rowRanges)
rownames(GSM_mat) <- paste0(tmp_name$seqnames , ":" , tmp_name$start , "-" , tmp_name$end)

## 转化为数值矩阵
GSM_mat_num <- apply(GSM_mat , 1 , as.numeric)
rownames(GSM_mat_num) <- colnames(GSM_mat)
GSM_mat_num <- t(GSM_mat_num)

## motfi
motif_se <- getMatrixFromProject(atac_proj, useMatrix="MotifMatrix")
## The deviations are the bias corrected deviations in accessibility. For each motif or annotation (rows), 
## there is a value for each cell or sample (columns) representing how different the accessibility for peaks with that motif or annotation 
## is from the expectation based on all cells being equal, corrected for biases.
motif_mat <- assays(motif_se)$deviation
## 转化为数值矩阵
motif_mat_num <- apply(motif_mat , 1 , as.numeric)
rownames(motif_mat_num) <- colnames(motif_mat)
motif_mat_num <- t(motif_mat_num)

## 每个细胞分布算peak和motif的激活程度
result_peak <- c()
for( cell in colnames(GSM_mat_num) ){
    peak_num <- length(which(GSM_mat_num[,cell]>0))
    motif_num <- length(which(motif_mat_num[,cell]>0))
    tmp <- data.frame( cell = cell , peak_num = peak_num , motif_num = motif_num )
    result_peak <- rbind( result_peak , tmp )
}

result_peak$cell_type <- cell_type
result_peak$cluster <- cluster

out_file <- paste0(out_path, "/cluster" , cluster , ".motif_peak.tsv")
write.table( result_peak , out_file , row.names = F , sep = "\t" )